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ABSTRACT 


The variational method of undetermined multipliers is used to derive a 
multivariate model for objective analysis. The model is intended for the 
assimilation of three-dimensional fields of rawinsonde height, temperature and 
wind and mean level temperature observed by satellite into a dynamically 
consistent data set. Relative measurement errors are taken into account. The 
dynamic equations are the two nonlinear horizontal momentum equations, the 
hydrostatic equation, and an integrated continuity equation. The model 
Euler-Lagrange equations are eleven linear and/or nonlinear partial 
differential and/or algebraic equations. A cyclical solution sequence is 
described. 

Other model features include a nonlinear terrain-following vertical 
coordinate that eliminates truncation error in the pressure gradient terms of 
the horizontal momentum equations and easily accommodates satellite observed 
mean layer temperatures in the middle and upper troposphere. A projection of 
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the pressure gradient onto equivalent pressure surfaces removes most of the 
adverse impacts of the lower coordinate surface on the variational adjustment. 


An evaluation of the variational objective analysis model appears in the 
following companion paper. 

1. Introduction 

The accurate diagnosis the state of the atmosphere with data measured 
from space-based platforms requires that these data be merged with data 


collected 

routinely 

by traditional methods . 

Because of 

diversities 

in 

variables 

measured, 

measurement accuracy, 

observation 

location. 

and 


observation density, it becomes necessary to do the blending in a way that the 
final analyses is at least an approximate solution of some mathematical 
relationship between the variables. 

The creation of a multivariate data set that is consistent with the 
mathematical expressions for the behavior of air in the free atmosphere is 
more typically accomplished through some form of data assimilation coupled 
with an initialization for a numerical prediction model. Observational data is 
blended with model forecast fields through an interpolation technique 
(Cressman, 1959; Gandin, 1963; Schlatter; 1975) in which the latter, used as a 
first guess, is updated with the observations. These methods are multivariate; 
wind observations are used in the interpolative analysis of the height and 
temperature fields and vice versa. Then an initialization procedure such as 
dynamic initialization (e.g., Miyakoda and Moyer, 1968; Nitta and Hovermale, 
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1969) or normal mode initialization (e.g., Baer and Tribbia, 1977; Machenhaur, 
1977; Daley, 1981) brings the hybrid data set into consistency with a 
numerical model equation set. Highly sophisticated dynamically consistent data 
assimilation schemes such as those described by McPherson, et al, , (1979), 
Bengtsson, et al., (1982), Ghil, et al. , (1979), Temperton (1984), and many 
others produce accurate representations of the state of the synoptic scale 
atmosphere. These hybrid data sets, though modified to be consistent with the 
scales of motion permitted by the models, also may be useful for diagnostic 
studies . 

This paper describes a nonlinear multivariate objective analysis for 
which the constrained state is defined by primitive equations. The purpose of 
this study is to achieve the dynamically constrained merger of satellite data 
with traditional data and to assess the impact of the former upon the 
analysis. The method is diagnostic. It does, however, require estimates for 
the explicitly formulated local tendencies of the velocity components and 
temperature which might be provided by a numerical model. The completed 
diagnostic scheme will contain complex mathematical constraints that are 
similar to those found in sophisticated methods for data assimilation and 
initialization of numerical models. This study focuses on the diagnostics of 
cyclone systems, including precipitation systems within cyclone systems, 
whereas many studies of numerical forecasts with mixed data sets are focused 
upon improved forecast skill. However, the completed model may eventually be 
of value in comparisons with existing data assimilations. The basic model 
(MODEL I) derived through methods of variational calculus is described herein 
and a companion paper (Achtemeier, et al. , 1988) deals with an evaluation of 
the method using a case study. 
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The goal of our research is the development of a variational data 
assimilation method that incorporates as dynamical constraints, the primitive 
equations for a moist, convectively unstable atmosphere and the radiative 
transfer equation. Variables that can be included in the adjusted are the 
three-dimensional vector wind, height, temperature, and moisture from 
rawinsonde data, and cloud-wind vectors, moisture, and radiance from satellite 
data. This presents a formidable mathematical problem. In order to facilitate 
thorough analysis of each of the model components, we defined four variational 
models that divide the problem naturally according to increasing complexity. 
The first variational model (MODEL I) contains the two nonlinear horizontal 
momentum equations, the integrated continuity equation, and the hydrostatic 
equation. Problems associated with an internally consistent finite difference 
method, a nonlinear hybrid terrain-following vertical coordinate, formulations 
for the pressure gradient terms, formulations for the velocity tendency terms 
and the development of a convergent solution sequence are addressed with MODEL 
I and are the subject of this paper. 

MODEL II contains MODEL I plus the thermodynamic equation for a dry 
adiabatic atmosphere. The introduction of this additional constraint violates 
the requirement that the number of subsidiary conditions (dynamic constraints) 
must be at least one less than the number of dependent variables (Courant, 
1936). Inclusion of the same number of constraints as dependent variables 
overdetermines the problem and a solution is not guaranteed. Therefore, we 
will design MODEL I to increase the number of dependent variables. MODEL III 
contains MODEL II plus an additional moisture variable and an equation to 
describe moist adiabatic processes. MODEL IV includes MODEL III plus radiance 
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as a dependent variable and the radiative transfer equation as a constraint. 


The next section presents the philosophy of a variational objective 
analysis model. In Section 3 appear the dynamic equations in the forms they 
enter the variational formalism as constraints. The variational equations are 
derived in Section 4. Details concerning the grid mesh, boundary conditions, 
and convergence of the equations are found in Section 5. Section 6 summarizes 
the model. 

2. A Variational Approach to Diagnostic Data Assimilation 

Our diagnostic objective analysis is an adaptation of Sasaki 1 s (1958) 
method of variational analysis. Data from different measurement systems are 
weighted according to measurement accuracies and are blended using a least 
squares method into a hybrid data set that satisfies a set of subsidiary 
conditions. Sasaki (1970a) presented two variational formulations for the 
solution of the data assimilation problem. His "weak constraint" formalism 
requires only a partial satisfaction of the subsidiary conditions through 
coefficients determined by the analyst. The subsidiary conditions are 
satisfied exactly through the "strong constraint" method. Ikawa (1984) has 
shown that the weak constraint algorithm converges to the strong constraint 
formalism as the coefficients become large. 

This study makes use of the method of undetermined multipliers (strong 
constraint formalism). The constraints are the nonlinear horizontal momentum 
equations (products of a variational principle (Wang, 1984)), the hydrostatic 
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equation and an integrated form of the continuity equation. The adjustments 
are carried out on fields of meteorological variables obtained through 
univariate objective interpolation. This kind of variational formulation has 
been criticized by Williamson and Daley (1983) on the grounds that the 
adjustments to the dynamic state are carried out from gridded fields rather 
than from the observations. Alternatively, observation statistics for 
different measurements of the same variable can be carried in the analyzed 
fields, perhaps as proposed by Baker (1983). Another approach would blend the 
variational analysis with a optimized successive corrections methods similar 
to that proposed by Bratseth (1986). 

Another difficulty with the method of undetermined multipliers is the 
complexity of the variational equations which stimulates a need for simpler 
methods to create hybrid, dynamically balanced data sets. Wahba and 
Wendelberger (1980) have shown that multivariate statistical objective 
analysis and variational analysis are interchangeable for linear constraints. 
Our variational method permits nonlinear constraints, allows for the physical 
interpretation of the adjustments and provides mutual adjustment between the 
mass and wind fields (Temperton, 1984). 

Accurately gridded meteorological variables are a requirement for any 
good diagnostic analysis. There are also quantities which, because of poor 
instrument accuracy or insufficient sampling frequency, cannot be measured 
directly and must be inferred through functions of other measured variables; 
in our case, they are products of the variational blending process. Among 
these are hypersensitive variables such as vertical velocity and the local 
tendencies of the horizontal velocity components, that are sensitive to small 
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changes in the other variables. The variational diagnostic model must also 
produce accurate fields of these hypersensitive variables. 

Krishnamurti (1968) calculated diagnostic vertical velocities through a 
12-forcing function balance omega equation. More recently. Smith and Lin 
(1978) preferred vertical velocities diagnosed from the O'Brien (1970) 
variational method. Our variational model calculates vertical velocity from a 
generalized form of the kinematic method for which the O'Brien method can be 
shown to be a special case. 

The local tendency terms of the horizontal velocity components are 
particularly difficult to determine accurately because the coarse sampling 
frequency of operational data collection networks is not sufficient to resolve 
the frequency of meteorological disturbances. Local tendencies can be 
incorporated into the variational analysis by fixing them and assuming that 
the generated error will not appreciably contaminate the solution. But this 
ignores the fact that the tendency terms are of the same order of magnitude as 
the advection terms and that generated error undoubtedly will contaminate the 
solution, especially the error sensitive divergence calculations. Sasaki 
(1970b), Sasaki and Lewis (1970), and Lewis and Grayson (1972) have used a 
"time-wise localized" method which physically is not a time adjustment, but a 
space filter designed to adjust variables in space at a particular time such 
that the local tendency is minimized with partial constraint satisfaction. 
Achtemeier (1975) included local rates of change in a primitive equation 
variational model through a subsidiary variational formulation based upon 
O'Brien's (1970) divergence adjustment method. This method was considered a 
failure after an extensive analysis (Achtemeier, 1979) found unrealistically 
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large velocity component tendencies where actual velocity changes over a 12-hr 
period were small. 

More recently, Lewis (1980, 1982) has examined the problem of time 

consistency from a Lagrangian approach through the application of Thompson's 
(1969) variational method. Lewis et al. . (1983) and Lewis and Derber (1985) 

combined rawinsonde data with VAS height data taken 2.5 hr later and found 
vertical velocity fields that compared favorably with space-observed cloud 
fields and surface weather reports. These studies and the results from Bloom's 
(1983) mesoscale analysis imply that variational methods can be used with some 
success in the direct determination of tendency variables, at least for 
observation frequencies on the order of 3-hr. 

3. The Formulations for the Dynamic Constraints 

The dynamic constraints, for MODEL I are the two nonlinear horizontal 

momentum equations, the hydrostatic equation, and an integrated continuity 
equation. These constraints must undergo several transformations for the 
successful solution of the variational problem with the method of undetermined 
Lagrange multipliers. The equations are nondimens ionalized following the 
methodology of Charney (1948) and Haltiner (1971). These equations transform 
into a nonlinear terrain-following vertical coordinate. Another transformation 
removes a hydrostatic component introduced into lower coordinate surfaces by 
unlevel terrain and thermodynamic variables are partitioned to give reference 
and perturbation atmospheres. The reference atmosphere is dynamically 
consistent. The variational analyses are done on the perturbation variables. 


Following Shuman and Hovermale (1968), the horizontal momentum equation, 
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for the u-component (eastward directed component at a reference longitude) of 
the three-dimensional vector wind, nondimensionalized in an arbitrary vertical 
coordinate and Cartesian on a conformal projection of the earth is. 
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where conventional meteorological symbols are used. Anthes and Warner (1978) 
show that Y> a distortion that arises from the projection of points on a 
sphere on to a conformal surface, is at least two orders of magnitude smaller 
than f and can be neglected. The projection ratio, m, is defined for a Lambert 
conformal projection. As part of the nondimensionalization, the mapscale 
factor, m, and the Coriolis parameter, f, expand into m=l+R^K and f=l+R^C 
where the arrays K and C are of order one and R^=0.1. 


The nondimensionalized hydrostatic equation in the generalized vertical 
coordinate is , 


li . T 9 J- n P = o 

9a 3a 

a) Vertical Coordinate for the Variational Assimilation Model 


( 2 ) 


The vertical coordinate follows the concept of Phillips (1957). It blends 
from a terrain-following coordinate in the lower troposphere into a pressure 
coordinate in the middle troposphere. All horizontal variations with the lower 
coordinate surface are confined to levels below a reference pressure level p'\ 
This vertical coordinate has the following advantages for the variational 


assimilation model: 
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(1) The dynamical equations are presented in their simplest form on the 
pressure surfaces at and above p . Coding to omit terms that are zero for 
coordinate surfaces that are surfaces of constant pressure can result in a 
substantial reduction of computational overhead. The tradeoff is that the 
equations below p" are more complex than the equations written for the linear 
sigma coordinate. However, the magnitudes of these additional terms become 
small in the sigma levels above the lower coordinate surface. 

(2) All coordinate surfaces at and above p* are pressure surfaces. Vertical 
interpolation of initial fields of data, such as mean layer temperatures from 
satellite observations that function with pressure, is not required for these 
coordinate surfaces. Furthermore, there is no need to interpolate from sigma 
coordinates back to pressure coordinates in order to interpret the 
variationally adjusted fields of meteorological variables. 

(3) Where coordinate surfaces are not constant pressure surfaces, the pressure 
gradient terms of the horizontal momentum equations transform into two large 
and compensating terms where there is steep, sloping terrain. The variational 
formalism will separate the pressure gradient terms and combine the large 
uncompensated terms with terms from the other equations. The large 
nonmeteorological impacts by these terms can cause significant error in the 
final solution. The nonlinear vertical coordinate eliminates this problem for 
the middle and upper troposphere (coordinate surfaces above the elevation of 
p\) A partition to reduce the impacts by these terms for layers below p' v is 
the subject of section b). The smooth transition from the terrain-f ollowing to 
the pressure coordinate is accomplished by fitting two curves which are 
piecewise continuous through the second derivatives. The curve for the upper 
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layer bounded by p u at the top and by p" at the bottom is given by a straight 
line subject to the boundary conditions that o=0 at p=p u and that a=a" at 
p=p /x . This equation is 


a 


a* 


p-p 

r r u 
P*-P 

u 


(3) 


The equation for the nonlinear part of the hybrid vertical coordinate between 


p" and the surface 

pressure p s 

is 

found subject 

to the following four 

conditions: 





a = 

1.0 


at p=p s 


a = 

a* 

| 



do 

Sp 

1 

a*/(p*-p u ) 


at p=p* 

(4) 

3p" 





These four conditions specify the equation as a cubic polynomial which takes 


the form 
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6 (p-p*) 3 


+ a* 


( p-p u } 

(p*-p u ) 


(5) 
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Fig. 1 shows the relationship between sigma and pressure for the levels 
below the elevation of the 600 mb pressure surface for the coordinate 
parameters that have been selected for the variational assimilation. The 
reference pressure p" is at 700 mb. A straight line from 700 mb to 1000 mb 
separates two sets of curves which describe the relationship between sigma and 
pressure for low surface pressure (high elevation) from those for high surface 
pressure. Wherever the slopes of the curves in Fig. 1 are less than the slope 
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of the straight line, the pressure thicknesses between sigma coordinate 

surfaces are compacted. The smallest pressure thicknesses are found nearest 

the lower coordinate surface. These layer depths increase to approach the 

*•* 

thicknesses of the pressure layers at levels above p . 

The slopes of the curves in Fig. 1 exceed the slope of the straight line 
at locations where the surface pressure is greater than 1000 mb. The pressure 
thicknesses between the sigma coordinate surfaces have been expanded. Note how 
the curve from 1100 mb converges asymptotically into the straight line by 
p=900mb. The nonlinear vertical coordinate forces most of the transition 
between terrain-following coordinate surfaces and pressure-following 
coordinate surfaces into the layer immediately above the ground. Thus 
coordinate surfaces in the lower troposphere tend to behave as pressure 
surfaces that are punctuated by areas of higher elevation. 

Fig. 2 shows the distribution of coordinate surfaces below 600 mb as the 
surface pressure varies from 800 to 1025 mb, the approximate range of surface 
pressures for the smoothed orography of the variational model. The greater 
compression of the coordinate surfaces over higher elevation nearest the 
surface is clearly evident. Notice how the nonlinear coordinate surfaces tend 
to become surfaces of constant pressure at locations away from the areas of 
high elevation. Note also the increased pressure depth of the lowest layer 
where the surface pressures exceed 1000 mb. Clearly this nonlinear vertical 
coordinate does not provide for a boundary layer of uniform thickness. 
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b) Partition to Reduce Impact of Nonmeteorological Hydrostatic Terms 

Below p'' where coordinate surfaces are not constant pressure surfaces, 
the pressure gradient terms of the horizontal momentum equations transform 
into two large and compensating terms where there is steep, sloping terrain. 
The variational formalism separates these pressure gradient terms and combines 
the large uncompensated terms with terms from the other equations. The impacts 
by these terms can introduce significant error into the final solution. This 
problem can be avoided by partitioning the hydrostatic equation into terrain 
terms and meteorological perturbation terms and subjecting the latter to the 
variational operation. Note that what is done is a partitioning, not a 
transformation. There is no change in the vertical coordinate. 

Consider the geopotential height and temperature in (2) to contain the 
"whole" signal and partition them into terrain, reference, initial, and 
adjusted variables according to, 

( ) w - ( ) T +( ) R +( )° +[()-( )°] (7) 

Furthermore, partition the whole pressure into 

Pw = PT + Pe ( 8 ) 

where the subscript, e, identifies "equivalent" pressure surfaces as a 
distinction between the adjustable heights obtained by this method and the 
heights of the pressure surfaces from the original observations. Upon 
partitioning, the hydrostatic equation can be expressed as the sum of four 
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groups of terms. These are: 


Terrain, 


P 30 ° 3cJ>_ T 0 d In p^ 
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Adjustment, 


[(^ - 


1 r\ 


(ID 


Perturbation, 


r 30 

[9J + y T ] (12) 

where Y =kln(p e )/^(T . Since the subscript w terms are known from the 
observations, the terrain height and its vertical gradient may be found from 
the specification of p e and requiring the terrain group to sum to zero. An 
accurate solution depends upon specification of a representative pressure for 
each layer. After some experimentation, it was found that, if the equivalent 
pressure at the top and the bottom of the layer is known, the average of the 
arithmetic mean plus twice the geometric mean, 

p - 3 tj ( p 2 + p x ) + 2 (13) 

gives an accurate representative pressure. 


Fig. 3a shows the heights of the lower coordinate surface at 1200 GMT 10 
April 1979 without partitioning of the hydrostatic equation. The figure shows 


mostly the variation of elevation. Fig. 3b shows the meteorological heights 
analyzed to the 1000 mb pressure surface. The low^ about 100 mb deep over 
Colorado^ is an order of magnitude smaller than the height of the smoothed 
terrain. Fig. 3c shows the heights of the 1000 mb equivalent pressure surface 
after partitioning the hydrostatic equation. The resemblance of all features 
to the heights of the actual 1000 mb surface is evident except for the higher 
central height of the low center over Colorado. The underestimation of this 
feature occurred because colder temperatures at higher elevations (the terrain 
temperature) were not partitioned. Since we have merely partitioned the 
heights, not neglected height terms, the remaining heights that make up the 
difference in the heights between Fig. 3c and Fig. 3b are spread among the 
remaining three groups of terms that make up the hydrostatic equation. 

Upon removal of the terrain height, the heights are averaged over each 
level to obtain the reference height. Then a hydrostatic reference atmosphere 
is found by requiring the reference group of terms to sum to zero. 

The residual group includes the small part of the terrain group that is 
subject to variation. These terms tend to compensate. Since the variational 
adjustments are only few tenths to one degree Kelvin, the residual terms are 
at least two orders of magnitude smaller than the comparable terrain terms and 
are one order of magnitude smaller than the perturbation terms. If the 
adjustment terms are represented by 3> then the hydrostatic equation is given 
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Once the hydrostatic equation is partitioned, the pressure gradient terms 
of the horizontal momentum equations can be partitioned to reduce the impacts 
of unlevel terrain. The partitioned pressure gradient term from (1) is, 


PGX 
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where , 
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is the sum of a small 
terrain terms, 
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c) Partition of the Local Tendencies 


Local changes in the horizontal velocity components are caused by a 
combination of translation of existing disturbances and development. In 
partitioning the tendencies, we note that, for example, the local change in 
the u- component of the wind caused by a moving weather system is 


3u 

3t 


• Vu 


+ 

dt 
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where c is the velocity of an advective or steering current (Fjortoft , 1952) , 
usually a smoothed middle tropospheric wind. Let u=uq+u* where ug is the 
u-component of the steady state part of the circulation and u 1 is the 
u-component arising from development. Then, 
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3u 

at 


v “o + <5T - c ' Vu,) 


(19) 


The first term of (19) is the local change in u caused by translation of the 
steady state part of the weather disturbance. The second term contains the 
local change of u from development. Note that the vertical advection of u is 
considered part of development. 


The use of the advective current throughout the troposphere is valid 
because most synoptic systems tend to maintain vertical structure. Any changes 
in vertical structure are assumed to be the result of development. The 
variational formalism requires that the adjustments be carried out on the 
total velocity components. Therefore, we represent the local tendency of u by 
(18). The total derivative, an approximate developmental component, is defined 
as a new dependent variable, e u =du/dt (e v =dv/dt). 

d) Transformation of the Integrated Continuity Equation 


The mass continuity equation in generalized coordinates (Shuman and 
Hovermale, 1968) is 
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The material derivative in the Lambert map projection is 
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Upon expanding the map scale factor m, (20) becomes 

„ ^ „ tr 2 ,9u/K , 9 v/K n , d /fl 9p x n 
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The last term of (22) is determined from the nonlinear vertical coordinate 
(5). Further, given (6) and the following definitions, 
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J s is obtained by substituting p s for p in (24). Including these modifications 
leads to the following form for the continuity equation, 


3u dv do • 
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where F = q 0) i + R K 2 + -^^) 

2 s 1 8x 3y ' 
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The q^ and q 2 arise through the nonlinearity of the coordinate transformation. 
They vanish where p>p . 


When solved for the vertical velocity, (28) becomes, 


a = a Q - HQ da 
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where 
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The integral function, Q, is proportional to the pressure thickness of the 
sigma layers. In order to simplify the Euler-Lagrange equations derived in the 
following section, we require Q=l. This assumption removes the dependence of 
the integrated divergence on the variable pressure thickness of the sigma 
layers. Therefore, divergences in the layers near the surface over elevated 
terrain receive proportionally greater weight in the vertical velocity 
adjustment. 


e) Finite Difference Equations for the Dynamic Constraints 

Two variational models were derived with finite difference versions of 
the dynamic constraints. The first was derived with the dynamic equations 
written in uncentered differences on a uniform grid and the second was 
formulated with centered differences on a staggered grid. We sought difference 
formulations for the Euler-Lagrange equations that were symmetric about the 
central grid point. The centered difference formulation on a staggered grid 
proved most suitable from this standpoint. 


Shuman and Hovermale (1968) and Anthes and Warner (1978) define the 
horizontal finite difference operators and the finite averaging operators as 

a x = (0t i+l/2,j “ a i-l/2,j ) 1 Ax 
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The i is the east-west index, the j is the north-south index as measured at 
the grid origin which is located at the lower left corner of the grid. In 
addition, the vertical differences and averages are defined by 

% = <Vl/2 ' Vl/2> ' A0 

“ = <C W/2 + “k-1/2* 11 

Figure 4 shows the staggered grid developed for this model. The 
geopotential $ is defined at the grid intersections, v is located at the top 
and bottom and u is located at the sides of the grid square. The divergence D 
is found at the center of the grid. The layer mean temperatures T are defined 
at one half grid length above and below the grid intersections and the 
vertical velocity^ 6^ is located one half grid space above and below the 
divergence. Mesinger and Arakawa (1976) have shown that phase speed and 
dispersion properties of this staggered grid make it inferior relative to 
other grid configurations for numerical prediction. However, the grid with v 
located on the top and bottom and u located on the sides of the grid box is 
well suited for the solution sequence developed for the Euler-Lagrange 
equations in the following section. Other variables used in the variational 
analysis are collocated with the variables in Fig. 4 as follows: e v and A, at 
v , e u and ^ at u, ^at D, and at T . 


The finite difference equations for the horizontal momentum equations 
written for the staggered grid are 
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The analogues for the continuity and hydrostatic equations are 
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The four dynamic constraints are referenced, respectively, to the following 
locations: at v, M 2 at u, M 3 at D, and M A at T on Fig. 4. 


4. Euler-Lagrange Equations 


The variational analysis melds satellite data with conventional data at 
the second stage of a two-stage objective analysis. All data are gridded 
independently in the first stage and are combined in the second stage. Future 
versions will combine the two stages. The gridded observations to be modified 
are meshed with the dynamic constraints through Sasaki's (1970a) variational 
formulation. The finite difference analog of the adjustment functional is 
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The weights tt^, i=l,7 are Gauss* precision moduli (Whittaker and Robinson, 
1926). The gridded observations (u°, v°, 6°, <fc°, T°, e u °, e v °) to be adjusted 
enter in a least squares formulation and receive precision modulus weights 
according to their relative observation accuracies. The strong constraints to 
be satisfied exactly are introduced through the Lagrangian multipliers 


i=l,4. 


Objectively modified meteorological variables are determined by requiring 
the first variation on F to vanish. A necessary condition for the existence of 
a stationary set is that the functions are determined from the domain of 
admissible functions as solutions of the Euler-Lagrange equations. The 
variation is to be carried out at every point (r,s) within the grid. Thus, 
setting the weights a^ = bj =1 and differentiating the integrand (40) with 
respect to the arbitrary variable ( a r>s ^» the Euler-Lagrange operator in 
finite differences is 
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The delta functions 6*, 6g, equal 1 where r=i or s=j and are zero 

x 

elsewhere. Each term in ^i,j that contains an overbar term, e.g. a r>s > 
produces a corresponding overbar term in the Euler-Lagrange equations when 
subjected to the operations specified by (41). It is convenient that the 
multiplicate overbar terms such as a r>s that appear in the nonlinear terms of 
the constraints be replaced by a r>s so that fewer gridpoints are required to 
express these terms in the Euler-Lagrange equations. 
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The Euler-Lagrange equations for u, v, and a result from the operations 
specified by (41). The equations are 
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The Euler-Lagrange equations for the thermodynamic variables $ and T are 
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Similarly, the operations performed for c u and e v yield 
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Variation on the Lagrange multipliers restores the four original dynamic 
constraints (35)-(38). 

Some of these Euler-Lagrange equations are complicated nonlinear partial 
differential equations for which solutions are difficult to obtain by direct 
methods. We observe, however, that the nonlinear terms are all products with 
the Rossby number or with . These equations may therefore be linearized and 
a solution obtained through a cyclical method as follows. Terms multiplied by 
R Q or R-l are expressed with observed variables at the first cycle, and are 
expressed by previously adjusted variables at higher cycles. At any particular 
solution cycle, these terms and the terms that are determined by observed 
variables are specified and can be treated as forcing functions. We emphasize 
that this solution method is valid only for the latitudes and scales of motion 
where the Rossby number is less than one. 


Upon linearization, the Euler-Lagrange equations for u, v, a, and <I> 
become 
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Similarly, the four dynamic constraints become 


<t> - v + F r =0 

x 5 
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Now these equations and (46)-(48) complete a set of eleven simple 
algebraic or linear partial differential equations. Variables may be easily 
eliminated to reduce the number of equations. Equations (49), (50), (53), and 
(54) formulated as vorticity expressions are combined to eliminate u and v. 

A + A_ = (II <{> ) + (IL <f> ) — F, + F_ + (ILF.) + (ILF.) 

lx 2y 1 x x 1 y y ly 2x 1 5 x 1 6 x (.5/) 

Equation (56) is combined with (46) to eliminate T, 

IL 
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Note that both (57) and (58) contain terms that obey the identity 


(AB ) = A Z B + A B Z 

z z zz z z 


(60) 


Now , )) 2 » and can be eliminated in (52), (57), and (58). This leaves a 
three-dimensional second-order partial differential equation with non-constant 
coefficients in <t>; 
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where 


- F. + F„ + 
ly 2x 


< n i F 5>x + 


< n iVy 


+ F, + F, 


(62) 


All of the coefficients of the geopotential and its derivatives are functions 
of sums and/or derivatives of precision modulus weights. The coefficients of 
the three second order partial derivatives are sums of precision moduli. These 
are always positive because the precision moduli are always positive. Thus 
(61) must be everywhere elliptic. 


We also derive a diagnostic equation in A3. First, (49) and (50) are 
divided by and reformulated as components of the divergence. Then they are 
combined into the divergence and integrated in the vertical to yield. 
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Now a is eliminated from (51) and (55) and the result replaces the integrated 
divergence in (63). Application of the identity (60) to several terms of (63) 
yields the two-dimensional second-order elliptic partial differential equation 
for )> 3 , 
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where 


10 


^ < X 2 + Vh + ^ <- X l + W do ‘if + F 7 


The relationship between u, v, and A 3 in (49) and (50) shows A 3 is an 
adjustment velocity potential. Equation (64) must be solved for the adjustment 
velocity potential. With the exception of a few small terms that contain the 
divergent part of the wind, (61) is a diagnostic equation for the rotational 
part of the wind. All other dependent variables can be calculated once $ and 
^3 have been determined. The solution sequence requires that the divergent 
components enter the velocity components, calculated from (53) and (54), 
through forcing functions that depend upon "old" data. In order for the 
divergent part of the wind to be current with the rotational part of the wind, 
the adjusted wind components are first found through (53) and (54) and these 
are readjusted for the divergent part of the wind. This solution method is not 
strictly rigorous, requiring an adjustment of an adjustment. Its verification 
through a case study is the subject of the accompanying paper. 


5. Computational Details 


The ten level variational assimilation model has the state variables 
staggered in both the horizontal and vertical dimensions. See Fig. 4 for the 
horizontal grid template. The variables u, v, e u , e v , and $ are located at 100 
mb intervals from the top of the domain (100 mb) to 700 mb (p”'). The 
constraints M 3 , M 2 , and M 3 are referenced to these surfaces. T, a, and M^ 
appear at 150- , 250- , 350- , 450- , 550- , and 650-mb surfaces. Further, the 
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upper boundary on a is at 50 mb (a=0). Below p' c , u, v and the developmental 
components appear on sigma surfaces and a and are located at the half 
levels. The first three dynamic constraints and $ are referenced to the 
equivalent pressure levels and the mean layer temperature is located at the 
half levels. The lower boundary for 6 (6=0) is the ground. We have also chosen 
the surface observations to be representative of the average conditions of the 
lowest sigma layer. This means that the boundary layer divergence is 
representative of the mean divergence of this lowest layer. 

The correct number of boundary conditions are furnished by the 
variational formulation such that a unique solution is provided when natural 
and/or imposed boundary equations are satisfied (Forsythe and Wasow, 1960). 
Natural boundary conditions are derived from the constraints as numerical 
expressions to be solved. However, because the MODEL I dynamic constraints 
produce complex boundary conditions, we use imposed boundary conditions and 
specify the dependent variables on the boundaries. The requirements of the 
solution method for MODEL I are that boundary conditions be specified for the 
geopotential and adjustment velocity potential diagnostic equations, the 
remaining equations in the adjustment cycle, and the vertical velocity. In 
addition, • special boundary conditions are imposed by the cyclic solution 
sequence. Details of the various boundary conditions follow. 

a) Boundary Conditions on the Geopotential Adjustment 

Imposed boundary conditions for the geopotential adjustment equation (61) 
are supplied by the gridded fields of the observed meteorological component of 
the partitioned height field. The top boundary is provided by the analysis at 
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100 mb. The lower boundary is the meteorological height component transformed 
onto equivalent pressure surfaces. 

b) Boundary Conditions on the Velocity Adjustment Potential 

Lateral boundary conditions are required for the adjustment velocity 
potential equation (64). It is well known that the specification of boundary 
conditions on the velocity potential determines the structural details of the 
recovered wind field to some degree (Hawkins and Rosenthal, 1965). 
Furthermore, there appears to be no method of uniquely specifying the boundary 
conditions (Shukla and Saha, 1974; Eskridge, 1977; Liu, 1977, Stephens and 
Johnson, 1978). Dirichlet boundary conditions force all of the divergence 
adjustment into the v-component along the x-boundaries and into the 

u- component along the y-boundaries. Neumann boundary conditions force the 

adjustment into the u- component along the x-boundaries and into the 

v-component along the y-boundaries. We used Dirichlet conditions and set 
the adjustment velocity potential equal to zero on all boundaries. This choice 
for boundary conditions left the divergent part of the wind unadjusted at the 
boundaries . 

c) Boundary Conditions on the Vertical Velocity 

The boundary conditions on a are 6- 0 at the ground and at 50 mb. Because 
the lower coordinate surface slopes with the underlying terrain, there may 
exist vertical velocity near the ground. Given as 0> s =dp s /dt, the surface 
vertical velocity is a combination of flow over elevated terrain and through 
evolving meteorological pressure fields. Upon partitioning into terrain and 
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meteorological components, the surface vertical velocity is 

9p m 

“s = - (V * 7p T + ^r + V * V P m ) (65) 

Scale analysis of the terms of (65) showed the first term in parentheses to be 
at least an order of magnitude larger than the meteorological terms. We 
therefore approximated the surface vertical velocity with the first term. 

d) Boundary Conditions on the Remaining Variables 

Horizontal boundary conditions are not required to determine the 
remaining variables in the interior of the analysis domain. However, boundary 
values of these variables are required in order to calculate horizontal 
derivatives for the forcing functions in subsequent solution cycles. Interior 
fields are extrapolated across the boundaries by using an approximation that 
is the sum of a locally averaged curvature with one half of a locally averaged 
gradient. This method provides boundaries that are compatible with the 
adjusted fields; they are generated, however, to eliminate boundary 
discontinuities, and do not satisfy the dynamic equations. 

e) Boundary Conditions Required by the Cyclic Solution Method 

Initial tests of the variational assimilation model with the case study 
described in the following article (Achtemeier, et al. , 1988) revealed local 

violations of linear stability along the lateral boundaries. These 
instabilities spread into the interior of the model domain with successive 
cycles. The adjustment of the geopotential height field (61) is forced to take 
on the gridded values of the observed geopotential at the boundaries 
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regardless of the relative weights ascribed to the other variables. Small 
perturbations in the heights near the boundaries are frozen into the 
geopotential adjustment and cause the instabilities. We were unable to totally 
eliminate the perturbations but were able to eliminate the buildup of the 
undesired waves by requiring variational analysis to satisfy the geostrophic, 
hydrostatic equations near the boundaries. These solutions grade into the 
solutions for the full nonlinear dynamic equations at five grid spaces into 
the grid interior. 


f) Convergence Criteria 


The convergence criteria for the general second-order partial 
differential equation with nonconstant coefficients, 

aX x» + bX yy + cX OT + dX x + eX y + - gX + h ‘ ° (66) 

obtained by the partial wave technique is 
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The coefficients a, b, c, and g in the geopotential adjustment equation are 
always positive. Further, the coefficient d is just the horizontal derivative 
of a, e is the horizontal derivative of b, and f is the vertical derivative of 
c. The most stringent convergence requirement is that the absolute magnitude 
of the derivative of a coefficient not exceed the value of the coefficient. 
This requirement can be easily satisfied through the definitions of the 
precision modulus weights. The coefficients of the adjustment velocity 
potential are similarly related except that c=f=0. 
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Convergence of the cyclic solution sequence for MODEL I is assured for 
the latitudes and scales of motion where the Rossby number is less than one. 
When the Rossby number exceeds one, the adjustment terms in the forcing 
functions approach or exceed the magnitudes of the variables being solved for, 
a condition that favors the development of linear instability. A determination 
of the limits of convergence MODEL I is the subject of continuing research. 

6. Some Concluding Remarks 

We have presented an outline of the first of four variational 
assimilation models that will meld data collected from rawinsonde (wind, 
temperature, height, moisture) with data collected from space-based platforms 
(cloud wind vectors, moisture, mean-layer temperatures). This method is, by 
design, independent of numerical prediction models. The first model, MODEL I, 
incorporates as dynamical constraints, the two nonlinear horizontal momentum 
equations, the hydrostatic equation, and an integrated continuity equation. 
The vertical coordinate minimizes the interpolation from pressure to 
terrain-following coordinates, more easily accommodating mean-layer 
temperature data observed by satellite in the middle and upper troposphere, 
and decreases truncation error associated with the pressure gradient force in 
the horizontal momentum equations. Reformulations for the horizontal 
tendencies of u and v are designed to increase the accuracy of the variational 
analysis for these hypersensitive quantities. 

We designed a cyclical solution sequence that linearizes the eleven 
Euler-Lagrange equations that comprise MODEL I as functions of the Rossby 
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number. This solution method does not exclude ageostrophic motion, it only 
requires that the Rossby number be less than one. A potential problem in the 
calculation of the adjusted velocity components is that this solution method 
is not strictly rigorous, requiring an adjustment of an adjustment. Its 
verification through a case study is the subject of the following paper 
(Achtemeier et al. , 1988). 
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Figure Captions 


Figure 1. The relationship between sigma and pressure for the levels below the 
elevation of the 600 mb pressure surface for the coordinate parameters 
that have been selected for the variational objective analysis. 

Figure 2. The distribution of hybrid coordinate surfaces below 600 mb as the 
surface pressure varies from 800 to 1025 mb. 

Figure 3. Heights at the lower coordinate surface for a) nonpartitioned 

terrain-following coordinate, b) the 1000 mb pressure surface, and c) the 
equivalent pressure surface. 

Figure 4. A portion of the staggered grid used for the variational diagnostic 
model. 







